home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-20 / usat92b2.zip / USATQB.BAS < prev    next >
BASIC Source File  |  1993-02-03  |  28KB  |  890 lines
  1. ' AMSAT ORBITAL PREDICTION PROGRAM de W3IWI - May 1980
  2. ' COPYRIGHT 1980 by Dr. Thomas A. Clark, W3IWI
  3. '                   6388 Guilford Road
  4. '                   Clarksville, MD. 21029
  5.  
  6. ' REVISED & MODIFIED FOR IBM-PC by R. D. Welch, W0SL - May 1983
  7. '                                  908 Dutch Mill Dr.
  8. '                                  Manchester, Mo. 63011
  9.  
  10. ' ENHANCED AND DEBUGGED BY Ing. H.F.STRASSER, OE1HSI - JAN. 1985
  11. '                          A 1238
  12. '                          VIENNA, AUSTRIA
  13.  
  14. ' Modified to use the NASA/NORAD two-line element set format (as distributed
  15. ' by Dr. Kelso of the Air Force Institute of Technology) and revised for
  16. ' QBASIC by Antonio Querubin, AH6BW
  17. '    Internet:  tony@mpg.phys.hawaii.edu
  18. '    AMPRnet:  ah6bw@uhm.ampr.org
  19. '    PBBS:  ah6bw@ah6bw.hi.usa.oc
  20. '    BITnet:  querubin@uhunix
  21.  
  22. ' Last update:  January 15, 1993
  23.  
  24. ' Permission granted for non-commercial use providing
  25. ' credit is given to the author(s), AMSAT and ORBIT Magazine.
  26.  
  27. DEFDBL A-Z
  28.  
  29. DECLARE FUNCTION gets$ ()
  30. DECLARE SUB getkey (k$)
  31. DECLARE SUB timezone (d0$, t0$, d1$, t1$, tzcor!)
  32. DECLARE FUNCTION juliantodate$ (day%, year%)
  33. DECLARE FUNCTION datetojulian% (d$)
  34. DECLARE FUNCTION trim$ (i$)
  35. DECLARE FUNCTION readtle% (tlefile$, name$(), epochyear%(), epochjulianday#(), inclination#(), epochRAAN#(), eccentricity#(), epochArgOfPerigee#(), epochMA#(), MeanMotion#(), epochrevolution&())
  36. DECLARE SUB readgrnd ()
  37. DECLARE SUB plotloc (longitude!, latitude!, xcoord%, ycoord%)
  38. DECLARE FUNCTION arctan# (y#, x#)
  39.  
  40. DEF fnyear% = VAL(RIGHT$(d$, 2))
  41. DEF fnhour% = VAL(LEFT$(t$, 2))
  42. DEF fnmonth% = VAL(LEFT$(d$, 2))
  43. DEF fnday% = VAL(MID$(d$, 4, 2))
  44. DEF fnminute% = VAL(MID$(t$, 4, 2))
  45. DEF fnsecond% = VAL(RIGHT$(t$, 2))
  46. DEF fnt$ (d%) = MID$(STR$(d% + 100), 3)
  47. DEF fnleap% (year%) = ((year% MOD 4) = 0)
  48. DEF fnradians (degrees) = degrees / 180# * pi#
  49. DEF fndegrees (radians) = radians / pi# * 180#
  50. DEF fnarcsin (sine) = ATN(sine / SQR(1 - sine * sine))
  51. ' Greenwich Mean Sidereal Time
  52. DEF fngmst# (YR%) = (99.6367# - .2387# * (YR% - 1989) + .9856# * INT((YR% - 1989) / 4#)) / 360!
  53.  
  54. CLEAR
  55.  
  56. maxsats% = 17
  57. DIM name$(maxsats%)
  58. DIM epochyear%(maxsats%), epochrevolution&(maxsats%), satl(maxsats%, 2)
  59. DIM epochjulianday#(maxsats%), inclination#(maxsats%), epochRAAN#(maxsats%)
  60. DIM eccentricity#(maxsats%), epochArgOfPerigee#(maxsats%)
  61. DIM epochMA#(maxsats%), MeanMotion#(maxsats%)
  62. DIM BeaconFrequency&(maxsats%), SemiMajor#(maxsats%)
  63. DIM sat%(7), pkt%(7), kep%(7) ' Image arrays (9 x 5 pixels)
  64.  
  65. DIM cc#(3, 2) ' Translation matrix from the satellite orbital plane coordinate
  66.              ' system to the celestial coordinate system (origin at geocenter)
  67.    
  68. ' ****** NUMERIC CONSTANTS ******
  69. pi# = 3.141592653589793# ' Value of PI
  70. r0# = 6378.14#           ' Earth's radius
  71. f# = 1# / 298.16#        ' 1/Earth flattening coef.
  72. g0# = 75369793000000#    ' GM of Earth in (orbits/day)^2/km^3
  73. g1# = 1.0027379093#      ' Solar/Sidereal time ratio
  74. c# = 299792.458#         ' Speed of light (km/sec).
  75.  
  76. ' String constants
  77. tlefile$ = "USAT.TLE"   ' Default NASA/NORAD element file
  78.  
  79. ' YOU WILL NEED THE FILES USATCGA.MAP, USAT.TLE AND USATGRND.DAT
  80. ' TO RUN THIS PROGRAM
  81. ' **** SATAUS Menuprogramm de OE1HSI VERSION 1.5  26.JAN. 1985
  82.  
  83. ' ****** ESTABLISH SATELLITE ELEMENT MATRIX ******
  84. PRINT "Reading satellite elements from "; tlefile$
  85. numsats% = readtle%(tlefile$, name$(), epochyear%(), epochjulianday#(), inclination#(), epochRAAN#(), eccentricity#(), epochArgOfPerigee#(), epochMA#(), MeanMotion#(), epochrevolution&())
  86.  
  87. ' Initialize beacon matrix and check age of satellite elements.
  88. FOR i% = 1 TO numsats%
  89.   BeaconFrequency&(i%) = 100
  90.   IF epochyear%(i%) < fnyear% - 1 THEN
  91.     PRINT "Warning:  Elements for satellite "; trim(name$(i%)); " are very old."
  92.     PRINT "You should update your "; tlefile$; " file."
  93.   END IF
  94. NEXT
  95.  
  96. ' Initialize the ground station data
  97. CALL readgrnd
  98.  
  99. ' ****** DRAW AND STORE SATELLITE SYMBOLS ******
  100. CLS
  101. SCREEN 2
  102. LINE (4, 0)-(4, 4)
  103. LINE (0, 2)-(8, 2)
  104. GET (0, 0)-(8, 4), sat%
  105. PUT (0, 0), sat%
  106. CLS
  107. LINE (4, 1)-(4, 3)
  108. LINE (3, 2)-(5, 2)
  109. GET (0, 0)-(8, 4), pkt%
  110. PUT (0, 0), pkt%
  111.  
  112. ' ***** Initializations done. *****
  113.  
  114. 50
  115. ' ****** DERIVED CONSTANTS ******
  116. ' These values depend on the ground station data.
  117. ' C$=GROUND STATION CALL+LOCATION STRING
  118. c$ = trim$(gsid$) + " " + trim$(gsloc$)
  119. pi2# = 2 * pi#
  120. SinGSLat# = SIN(fnradians(gslat!))    ' sin/cos of ground station latitude
  121. CosGSLat# = COS(fnradians(gslat!))
  122. SinGSLong# = SIN(fnradians(-gslong!))  ' sin/cos of ground station longitude
  123. CosGSLong# = COS(fnradians(gslong!))
  124. R9 = r0# * (1 - (f# / 2) + (f# / 2) * COS(2 * fnradians(gslat!))) + gsheight% / 1000
  125. L8 = ATN((1 - f#) ^ 2 * SinGSLat# / CosGSLat#)
  126. GSzcoord = R9 * SIN(L8)
  127. GSxcoord = R9 * COS(L8) * CosGSLong#
  128. GSycoord = R9 * COS(L8) * SinGSLong#
  129.  
  130. 80
  131. CLS
  132. SCREEN 0
  133. KEY(9) OFF
  134. KEY(10) OFF
  135. PRINT "                        USAT - QBASIC Version"
  136. PRINT
  137. PRINT "======================================================================"
  138. PRINT
  139. PRINT "            SELECT ONE OF THE FOLLOWING OPTIONS:"
  140. PRINT
  141. PRINT " (P) ORBITAL PREDICTION PROGRAM"
  142. PRINT
  143. PRINT " (R) REALTIME TRACKING AND HIGH RESOLUTION SCREEN"
  144. PRINT
  145. PRINT " (S) Display ELEMENTS OF SATELLITES"
  146. PRINT
  147. PRINT " (G) CHANGE OR ENTER GROUNDSTATION DATA"
  148. PRINT
  149. PRINT " (D) RETURN TO DOS"
  150. PRINT
  151. PRINT
  152. PRINT "ENTER SELECTION (P,R,C,G,D)--> "
  153. CALL getkey(k$)
  154. ON INSTR("PRSGD", k$) GOTO 4820, 2090, 240, 1620, 9000
  155. BEEP
  156. GOTO 50
  157.  
  158. 240
  159. ' ****** SATFILE.BAS - VERSION 1.0, ISSUE 1.0 - HSIMODIF.1/25/85 ******
  160. CLS
  161. DO
  162.   PRINT "Elements of the following satellites are available:"
  163.   PRINT
  164.   FOR j% = 1 TO numsats%
  165.     PRINT name$(j%),
  166.   NEXT
  167.   PRINT
  168.   ' ****** FIND SATELLITE RECORD ******
  169.   DO
  170.     PRINT : INPUT "Which satellite"; v1$: IF v1$ = "" THEN 80
  171.     v1$ = UCASE$(trim$(v1$))
  172.     FOR j% = 1 TO numsats%
  173.       IF UCASE$(trim$(name$(j%))) = v1$ THEN
  174.         CLS
  175.         ' ****** DISPLAY SATELLITE RECORD ******
  176.         PRINT "Satellite       = "; name$(j%)
  177.         PRINT "Epoch year      = "; epochyear%(j%)
  178.         PRINT "Epoch day       = "; epochjulianday#(j%)
  179.         PRINT "Inclination     = "; inclination#(j%)
  180.         PRINT "R.A.A.N.        = "; epochRAAN#(j%)
  181.         PRINT "Eccentricity    = "; eccentricity#(j%)
  182.         PRINT "Arg. of perigee = "; epochArgOfPerigee#(j%)
  183.         PRINT "Mean anomaly    = "; epochMA#(j%)
  184.         PRINT "Mean motion     = "; MeanMotion#(j%)
  185.         PRINT "Epoch orbit no. = "; epochrevolution&(j%)
  186.         PRINT "Beacon freq.    = "; BeaconFrequency&(j%)
  187.         PRINT
  188.         EXIT DO
  189.       END IF
  190.     NEXT
  191.     PRINT "Satellite not found."
  192.   LOOP
  193. LOOP
  194.  
  195. 1620
  196. ' ******* Groundsation data change v.1.0 OE1HSI  jan.-1985**********
  197. CLS
  198. PRINT "CURRENT GROUND STATION DATA"
  199. PRINT
  200. GOSUB 2080
  201. DO
  202.   PRINT "Do you want to CHANGE this DATA ? (Y/N)"
  203.   CALL getkey(k$)
  204.   ON INSTR("YN", k$) GOTO 1780, 2060
  205.   BEEP
  206. LOOP
  207. 1710
  208. PRINT
  209. PRINT
  210. GOSUB 2080
  211. DO
  212.   PRINT "Do you want to make further changes ? (Y/N) "
  213.   CALL getkey(k$)
  214.   ON INSTR("YN", k$) GOTO 1780, 2050
  215.   BEEP
  216. LOOP
  217. 1780
  218. PRINT
  219. PRINT "ENTER NEW DATA OR <RETURN> FOR UNCHANGED DATA":
  220. PRINT
  221. PRINT "Ground Station Identifier                    : ";
  222. LINE INPUT u$
  223. IF trim$(u$) <> "" THEN gsid$ = trim$(u$)
  224. PRINT "Location of station                          : ";
  225. LINE INPUT u$
  226. IF trim$(u$) <> "" THEN gsloc$ = trim$(u$)
  227. PRINT
  228. PRINT "Longitude WEST of Greenwich (max +360) or East of Greenw. entered as -0 to -180"
  229. PRINT
  230. INPUT "Enter (with decimals)                 : ", u$
  231. IF trim$(u$) <> "" THEN gslong! = VAL(u$)
  232. IF gslong! < 0 THEN gslong! = 360 + gslong!
  233. PRINT
  234. PRINT "LATITUDE NORTH of Equator + (max 90) SOUTH of Equator - (max 90)"
  235. PRINT
  236. INPUT "ENTER (With decimals)                 : ", u$
  237. IF trim$(u$) <> "" THEN gslat! = VAL(u$)
  238. PRINT
  239. INPUT "Station height above sea level in meters     : ", u$
  240. IF trim$(u$) <> "" THEN gsheight% = VAL(u$)
  241. PRINT
  242. PRINT "Correction from local computer time to UTC"
  243. PRINT
  244. INPUT "Enter (in hours with decimals)        : ", u$
  245. IF trim$(u$) <> "" THEN tzcor! = VAL(u$)
  246. GOTO 1710
  247. 2050
  248. PRINT
  249. PRINT "DATA SAVED AS USATGRND.DAT"
  250. OPEN "USATGRND.DAT" FOR OUTPUT AS #1
  251. PRINT #1, gsid$
  252. PRINT #1, gsloc$
  253. PRINT #1, gslong!, gslat!, gsheight%, tzcor!
  254. CLOSE #1
  255. GOTO 50
  256. 2060
  257. PRINT
  258. PRINT "DATA NOT CHANGED"
  259. GOTO 50
  260.  
  261. 2080
  262.   PRINT "Ground Station Identifier            : "; gsid$
  263.   PRINT "Location name                        : "; gsloc$
  264.   PRINT USING "West longitude (degrees)             : +###.##"; gslong!
  265.   PRINT USING "Latitude (degrees)                   :  +##.##"; gslat!
  266.   PRINT USING "Height above Mean Sea Level (meters) :   #####"; gsheight%
  267.   PRINT USING "Local time correction to UTC (hours) :  +##.##"; tzcor!
  268.   PRINT
  269. RETURN
  270. '**** END PROGRAM GROUNDSTATION DATA CHANGE/STORAGE OE1HSI  JAN. 1985 ****
  271.  
  272. 2090
  273. ' ****** ORBITS2 - VERSION 1.0, ISSUE 1.2 -11/1/83 ******
  274. KEY OFF
  275. SCREEN 2, 0
  276. CLS
  277.  
  278. ' Initialize D$ and T$ from date$ and time$
  279. CALL timezone(DATE$, TIME$, d$, t$, tzcor!)
  280. ' ****** SET UP KEY TRAPPING ******
  281. ON KEY(9) GOSUB 4760
  282. KEY(9) STOP
  283. ON KEY(10) GOSUB 4790
  284. KEY(10) OFF
  285. flg9 = 0
  286. flg10 = 0
  287. GOSUB 4290
  288.    
  289. ' ****** ORBIT DETERMINATION LOOP STARTS HERE ******
  290. q$ = ""
  291. DO UNTIL q$ = CHR$(27)
  292.   FOR j% = 1 TO numsats%
  293.     q$ = UCASE$(INKEY$)
  294.     IF q$ = CHR$(27) THEN EXIT FOR
  295.     CALL timezone(DATE$, TIME$, d$, t$, tzcor!)
  296.     t = INT(30.55 * (fnmonth% + 2)) - 2 * (INT(.1 * (fnmonth% + 7))) - 91
  297.     IF fnmonth% > 2 AND fnleap%(fnyear%) THEN t = t + 1
  298.     IF epochyear%(j%) = fnyear% - 1 THEN
  299.       t = t + 365
  300.       IF fnleap%(k%) THEN t = t + 1
  301.     END IF
  302.     t = t + fnday% + fnhour% / 24# + fnminute% / 1440# + fnsecond% / 86400#
  303.     GOSUB 6780
  304.     IF flg9 THEN CALL plotloc(longitude!, latitude!, plotx%, ploty%): GOSUB 3890
  305.     GOSUB 4040
  306.   NEXT
  307. LOOP
  308. GOTO 50
  309.    
  310. ' ****** PUT SATELLITE ON SCREEN ROUTINE ******
  311. 3890
  312.   GET (plotx% - 4, ploty% - 2)-(plotx% + 4, ploty% + 2), kep%
  313.   PUT (plotx% - 4, ploty% - 2), sat%, PRESET
  314.   PUT (plotx% - 4, ploty% - 2), sat%
  315.   PUT (plotx% - 4, ploty% - 2), kep%, PSET
  316.   PUT (plotx% - 4, ploty% - 2), sat%
  317.   IF FLG0 <> 0 THEN
  318.     PUT (satl(j%, 1), satl(j%, 2)), sat%
  319.     PUT (satl(j%, 1), satl(j%, 2)), pkt%, OR
  320.   END IF
  321.   satl(j%, 1) = plotx% - 4
  322.   satl(j%, 2) = ploty% - 2
  323.   IF j% = numsats% THEN FLG0 = 1
  324. RETURN
  325.  
  326. ' ****** PRINT SATELLITE DETAILS ROUTINE ******
  327. 4040
  328.   KEY(9) ON
  329.   KEY(10) ON
  330. 4050
  331.   KEY(10) STOP
  332.   KEY(9) STOP
  333.   IF FLGK = 1 THEN GOSUB 4270
  334.   IF flg9 = 0 THEN
  335.     LOCATE 2, 49
  336.     PRINT d$
  337.     LOCATE 2, 63
  338.     PRINT t$
  339.     IF elevation! >= 0 THEN
  340.       IF elevation! > 0 AND elevation! < 1 THEN BEEP
  341.     END IF
  342.     IF j% + 6 > 23 GOTO 4230
  343.     LOCATE j% + 6, 15
  344.     PRINT SPACE$(50)
  345.     LOCATE j% + 6, 8 - (LEN(trim$(name$(j%))) / 2)
  346.     PRINT trim$(name$(j%))
  347.     LOCATE j% + 6, 15
  348.   ELSE
  349.     IF flg10 = 1 THEN GOSUB 4540
  350.     IF flg10 <> 0 THEN
  351.       IF UCASE$(trim$(name$(j%))) <> u$ THEN RETURN
  352.       LOCATE 25, 69
  353.       PRINT "SELECTED";
  354.     END IF
  355.     LOCATE 25, 1
  356.     PRINT SPACE$(68);
  357.     LOCATE 25, 8 - (LEN(trim$(name$(j%))) / 2)
  358.     PRINT trim$(name$(j%));
  359.     LOCATE 25, 15
  360.   END IF
  361.   PRINT USING "###   ###  #####    #####  ###.#   ###.#    ######"; azimuth!; elevation!; range#; (r# - r0#); latitude!; longitude!; revolution&;
  362. 4230
  363.   IF flg9 <> 0 THEN LOCATE 20, 48: PRINT t$;
  364. RETURN
  365.    
  366. ' ****** SET UP SCREEN DISPLAY ROUTINE ******
  367. 4270
  368.   CLS
  369.   FLG0 = 0
  370.   FLGK = 0
  371. 4290
  372.   IF flg9 = 1 THEN
  373.     SCREEN 2, 0
  374.     DEF SEG = &HB800
  375.     BLOAD "USATCGA.MAP", 0
  376.     DEF SEG
  377.     CALL plotloc(gslong!, gslat!, plotx%, ploty%)
  378.     CIRCLE (plotx%, ploty%), 2
  379.     LOCATE , , 1, 12, 13
  380.     L1 = 22
  381.     L2 = 23
  382.     L3 = 24
  383.     LOCATE 20, 3
  384.     PRINT "Data for Groundstation             At Time=           UTC  On: "; d$
  385.     LOCATE 20, 26
  386.     PRINT gsid$
  387.   ELSE
  388.     ' FLG9=0
  389.     SCREEN 0, 1
  390.     LOCATE 1, 40 - LEN(c$) / 2, 0
  391.     PRINT c$
  392.     LOCATE 2, 5
  393.     PRINT "Real-time satellite tracking coordinates on            at          UTC"
  394.     LOCATE 25, 3
  395.     PRINT "F9 TOGGLES THE GRAPH/TABLE     F10 TO SELECT SINGLE SAT IN GRAPH   ESC TO END";
  396.     L1 = 4
  397.     L2 = 5
  398.     L3 = 6
  399.   END IF
  400.   LOCATE L1, 3
  401.   PRINT "  NAME OR   AZ    EL   RANGE   HEIGHT   LAT.   LONG.    ORBIT"
  402.   LOCATE L2, 3
  403.   PRINT "DESIGNATOR  DEG   DEG    KM      KM     DEG     DEG       NO."
  404.   LOCATE L3, 3
  405.   PRINT "----------  ---   ---  -----   ------  -----   -----    ------";
  406. RETURN
  407.    
  408. ' ****** SELECT SINGLE SATELLITE ROUTINE ******
  409. 4540
  410.   LOCATE 25, 1
  411.   PRINT SPACE$(79);
  412.   LOCATE 25, 1
  413.   INPUT ; "WHICH SATELLITE? (CR FOR ALL)"; i$
  414.   i$ = UCASE$(trim$(i$))
  415.   IF i$ = "" THEN
  416.     flg10 = 0
  417.   ELSE
  418.     u$ = i$
  419.     flg10 = 2
  420.   END IF
  421.   LOCATE 25, 1
  422.   PRINT SPACE$(79);
  423. RETURN
  424.  
  425. 4760 ' toggle FLG9
  426.   flg9 = -(flg9 = 0)
  427.   FLGK = 1
  428. RETURN 4050
  429.  
  430. 4790 flg10 = flg9
  431. RETURN 4050
  432.  
  433. 4820
  434.   '****** ORBIT2 - VERSION 2.0, ISSUE 1.0/HSI - 17/01/85 *****
  435.   KEY OFF
  436.   SCREEN 0
  437.   CLS
  438.   ' ****** HOUSEKEEPING ITEMS ******
  439.   C8$ = CHR$(10) + CHR$(10) + CHR$(10) + CHR$(10)
  440.  
  441.   DO
  442.     PRINT "Enter start date (mm-dd-yy):  ";
  443.     u$ = gets$
  444.     year% = VAL(RIGHT$(u$, 2))
  445.     month% = VAL(LEFT$(u$, 2))
  446.     day% = VAL(MID$(u$, 4, 2))
  447.     IF month% >= 1 AND month% <= 12 AND day% >= 1 AND day% <= 31 AND year% >= 91 AND MID$(u$, 3, 1) = "-" AND MID$(u$, 6, 1) = "-" THEN EXIT DO
  448.   LOOP
  449.   YY = year% + 1900
  450.   t$ = fnt$(year%) + "/" + fnt$(month%) + "/" + fnt$(day%) + " at "
  451.   D8 = day% + INT(30.55 * (month% + 2)) - 2 * (INT(.1 * (month% + 7))) - 91
  452.   IF month% > 2 THEN IF fnleap%(year%) THEN D8 = D8 + 1
  453.   PRINT "    Day #"; D8
  454.   PRINT
  455.   PRINT "Enter start time (hh:mm):  ";
  456.   DO
  457.     u$ = gets$
  458.     h% = VAL(LEFT$(u$, 2))
  459.     m% = VAL(RIGHT$(u$, 2))
  460.     IF h% >= 0 AND h% <= 23 AND m% >= 0 AND m% <= 59 AND MID$(u$, 3, 1) = ":" THEN EXIT DO
  461.     PRINT "Valid format is:  hh:mm"
  462.     PRINT "Re-enter:  ";
  463.   LOOP
  464.   T7 = D8 + h% / 24# + m% / 1440#
  465.   t$ = t$ + fnt$(h%) + fnt$(m%) + ":00 H"
  466.   PRINT "Enter duration in hours and minutes (hh:mm):  ";
  467.   DO
  468.     u$ = gets$
  469.     h% = VAL(LEFT$(u$, 2))
  470.     m% = VAL(RIGHT$(u$, 2))
  471.     IF h% >= 0 AND m% >= 0 AND m% <= 59 AND MID$(u$, 3, 1) = ":" THEN EXIT DO
  472.     PRINT "Valid format is:  hh:mm"
  473.     PRINT "Re-enter:  ";
  474.   LOOP
  475.   T8 = T7 + h% / 24# + m% / 1440#
  476.   PRINT USING "    From ###.####### to ###.#######"; T7; T8
  477.   DO
  478.     PRINT "Enter time increment (minutes):  ";
  479.     INPUT m%
  480.     IF m% > 0 THEN EXIT DO
  481.   LOOP
  482.   T9 = m% / 1440#
  483.   PRINT
  484.   INPUT "MINIMUM ELEVATION ? (DEFAULT 0) Deg. = ", E8
  485.   DO
  486.     PRINT
  487.     INPUT "Output to Printer (P) or Screen (S) ?-->", P$
  488.     P$ = UCASE$(P$)
  489.     IF (P$ = "P" OR P$ = "S") THEN
  490.       EXIT DO
  491.     ELSE BEEP
  492.     END IF
  493.   LOOP
  494.   IF P$ = "P" THEN outfile$ = "lpt1:" ELSE outfile$ = "scrn:"
  495.   CLOSE #2
  496.   OPEN outfile$ FOR OUTPUT AS #2
  497.   IF outfile$ = "lpt1:" THEN
  498.     PRINT "Put the printer on-line and align the page,"
  499.     PRINT "then press any key when ready."
  500.     DO
  501.       getkey (u$)
  502.       IF u$ <> "" THEN EXIT DO
  503.     LOOP
  504.   END IF
  505.  
  506. ' ****** SELECT SATELLITE FROM MENU ******
  507.   CLS
  508.   PRINT "SATELLITE SELECTION MENU"
  509.   PRINT
  510.   FOR j% = 1 TO numsats%
  511.     PRINT name$(j%),
  512.   NEXT
  513.   PRINT
  514.   DO
  515.     PRINT
  516.     INPUT "SELECT satellite >", y3$
  517.     y3$ = UCASE$(trim$(y3$))
  518.     FOR j% = 1 TO numsats%
  519.       IF UCASE$(trim$(name$(j%))) = y3$ THEN EXIT DO
  520.     NEXT
  521.     PRINT "Satellite not found."
  522.   LOOP
  523.   PRINT
  524.   PRINT "Doppler calculated for frequ. = "; BeaconFrequency&(j%); " MHz"
  525.   INPUT " Change frequency to:  (0 for default) ", d
  526.   IF d <> 0 THEN BeaconFrequency&(j%) = d
  527.   PRINT #2,
  528.   PRINT #2, "Orbital ELEMENTS for "; name$(j%)
  529.   PRINT #2,
  530.   PRINT #2, "Reference epoch = "; epochyear%(j%); " +"; epochjulianday#(j%)
  531.   PRINT #2, "Starting  epoch = "; year%; " +"; T7; " = "; t$
  532.   PRINT #2,
  533.   PRINT #2, "Parameter"; TAB(20); "Reference"; TAB(40); "Starting"
  534.   t = T7
  535.   IF epochyear%(j%) = year% - 1 THEN
  536.     t = t + 365
  537.     T8 = T8 + 365
  538.     IF fnleap%(k%) THEN t = t + 1: T8 = T8 + 1
  539.   END IF
  540.   ' Calculate start time parameters
  541.   GOSUB 6780
  542.   PRINT #2, "Orbit Number "; TAB(20); epochrevolution&(j%); TAB(40); revolution&
  543.   PRINT #2, "Mean Anomaly "; TAB(20); epochMA#(j%); TAB(40); fndegrees(MeanAnomaly#)
  544.   PRINT #2, "Inclination  "; TAB(20); inclination#(j%)
  545.   PRINT #2, "Eccentricity "; TAB(20); eccentricity#(j%)
  546.   PRINT #2, "Mean Motion  "; TAB(20); MeanMotion#(j%)
  547.   PRINT #2, "S.M.A.,km    "; TAB(20); SemiMajor#(j%)
  548.   PRINT #2, "Arg. Perigee "; TAB(20); epochArgOfPerigee#(j%); TAB(40); ArgOfPerigee#
  549.   PRINT #2, "R. A. A. N.  "; TAB(20); epochRAAN#(j%); TAB(40); RAAN#
  550.   PRINT #2, "Freq.,MHz    "; TAB(20); BeaconFrequency&(j%)
  551.   OldRevolution& = -1
  552.   OldJulianDay% = -1
  553.   
  554.   '****** COMPUTATION LOOP ******
  555. 6400
  556.   t = t + T9
  557.   julianday% = INT(t)
  558.   IF julianday% <> OldJulianDay% THEN
  559.     OldJulianDay% = -1
  560.     OldRevolution& = -1
  561.   END IF
  562.   GOSUB 6780
  563.   IF elevation! >= E8 THEN
  564.     IF julianday% <> OldJulianDay% OR revolution& <> OldRevolution& THEN
  565.       IF julianday% <> OldJulianDay% THEN
  566.         GOSUB 6690
  567.         OldJulianDay% = julianday%
  568.         PRINT #2, " U.T.C.    AZ    EL  DOPPLER   RANGE   HEIGHT  LAT.  LONG.  PHASE"
  569.         PRINT #2, "HHMM:SS   deg   deg     Hz       km      km    deg    deg   <256>"
  570.       END IF
  571.       PRINT #2, TAB(21); "- - - ORBIT #"; revolution&; "- - -"
  572.     END IF
  573.     OldRevolution& = revolution&
  574.     T4 = t - julianday%
  575.     s4 = INT(T4 * 86400)
  576.     H4 = INT(s4 / 3600 + .000001)
  577.     M4 = INT((s4 - H4 * 3600) / 60 + .000001)
  578.     s4 = s4 - 3600 * H4 - 60 * M4
  579.     doppler# = -BeaconFrequency&(j%) * 1000000 * Vr / c#
  580.     PRINT #2, fnt$(H4) + fnt$(M4) + ":" + fnt$(s4);
  581.     PRINT #2, USING "   ###   ###  #######"; azimuth!; elevation!; doppler#;
  582.     PRINT #2, USING "   #####   #####  ###.#  ###.#   ###"; range#; (r# - r0#); latitude!; longitude!; phase%
  583.   END IF
  584.   IF t < T8 THEN GOTO 6400
  585.   PRINT #2,
  586.   PRINT "Do YOU have another INQUIRY  ?  (Y/N) "
  587.   PRINT
  588.   PRINT "Else you return to the MAIN MENU !"
  589.   PRINT
  590.   DO
  591.     CALL getkey(k$)
  592.     ON INSTR("YN", k$) GOTO 4820, 6670
  593.     BEEP
  594.   LOOP
  595. 6670
  596.   CLOSE #2
  597. GOTO 80
  598.     
  599. '****** PAGE HEADER SUBROUTINE ******
  600. 6690
  601.   PRINT #2,
  602.   PRINT #2, c$; "   Lat.="; gslat!; "  W.Long.="; gslong!; "  Ht.="; gsheight%;
  603.   PRINT #2, TAB(14); "- - - Minimum Elevation = "; E8; "Deg. - - -"
  604.   PRINT #2,
  605.   PRINT #2, TAB(14); "- - - DAY #"; julianday%; "- - -  "; juliantodate$(julianday%, INT(YY)); "  - - -"
  606.   PRINT #2,
  607. RETURN
  608.  
  609. '****** ORBIT DETERMINATION AND UTILITY ROUTINES ******
  610.  
  611. ' Generate position data within the satellite orbital plane
  612. 6780
  613.   ' Calculate the Mean Anomaly (in radians) within the orbital plane
  614.   OrbitPos# = epochrevolution&(j%) + epochMA#(j%) / 360 + MeanMotion#(j%) * (t - epochjulianday#(j%))
  615.   revolution& = INT(OrbitPos#)
  616.   phase% = INT((OrbitPos# - revolution&) * 256)
  617.   MeanAnomaly# = (OrbitPos# - revolution&) * pi2#
  618.   ' Newton-Raphson calculation of Eccentric Anomaly
  619.   EA# = MeanAnomaly# + eccentricity#(j%) * SIN(MeanAnomaly#) + .5 * (eccentricity#(j%) * eccentricity#(j%)) * SIN(2 * MeanAnomaly#)
  620.   DO
  621.     SinEA# = SIN(EA#)
  622.     CosEA# = COS(EA#)
  623.     slope# = 1 - eccentricity#(j%) * CosEA#
  624.     Yerror# = EA# - eccentricity#(j%) * SinEA# - MeanAnomaly#
  625.     IF ABS(Yerror#) < .000001 THEN EXIT DO ELSE EA# = EA# - Yerror# / slope#
  626.   LOOP
  627.   ' Calculate the satellite position in the orbit plane
  628.   SemiMajor#(j%) = (g0# / (MeanMotion#(j%) * MeanMotion#(j%))) ^ (1 / 3)
  629.   ba2# = 1 - eccentricity#(j%) * eccentricity#(j%)
  630.   orbitX# = SemiMajor#(j%) * (CosEA# - eccentricity#(j%))
  631.   orbitY# = SemiMajor#(j%) * SQR(ba2#) * SinEA#
  632.   r# = SemiMajor#(j%) * slope#
  633.   ' Calculate relative orientation of the satellite orbit
  634.   PrecessFactor# = 9.95 * ((r0# / SemiMajor#(j%)) ^ 3.5) / (ba2# * ba2#)
  635.   SinInclination# = SIN(fnradians(inclination#(j%)))
  636.   CosInclination# = COS(fnradians(inclination#(j%)))
  637.   RAAN# = epochRAAN#(j%) - (t - epochjulianday#(j%)) * PrecessFactor# * CosInclination#
  638.   ArgOfPerigee# = epochArgOfPerigee#(j%) + (t - epochjulianday#(j%)) * PrecessFactor# * (2.5 * (CosInclination# * CosInclination#) - .5)
  639.   ' Generate the satellite orbital-plane to
  640.   ' celestial coordinates translation matrix
  641.   SinRAAN# = SIN(fnradians(RAAN#))
  642.   CosRAAN# = COS(fnradians(RAAN#))
  643.   SinAoP# = SIN(fnradians(ArgOfPerigee#))
  644.   CosAoP# = COS(fnradians(ArgOfPerigee#))
  645.   cc#(1, 1) = (CosAoP# * CosRAAN#) - (SinAoP# * SinRAAN# * CosInclination#)
  646.   cc#(1, 2) = -(SinAoP# * CosRAAN#) - (CosAoP# * SinRAAN# * CosInclination#)
  647.   cc#(2, 1) = (CosAoP# * SinRAAN#) + (SinAoP# * CosRAAN# * CosInclination#)
  648.   cc#(2, 2) = -(SinAoP# * SinRAAN#) + (CosAoP# * CosRAAN# * CosInclination#)
  649.   cc#(3, 1) = (SinAoP# * SinInclination#)
  650.   cc#(3, 2) = (CosAoP# * SinInclination#)
  651.   ' Calculate position in the celestial coordinate system (origin at geocenter)
  652.   celestialX# = orbitX# * cc#(1, 1) + orbitY# * cc#(1, 2)
  653.   celestialY# = orbitX# * cc#(2, 1) + orbitY# * cc#(2, 2)
  654.   celestialZ# = orbitX# * cc#(3, 1) + orbitY# * cc#(3, 2)
  655.   ' Now take Earth's rotation and solar orbit into account
  656.   GreenwichDeclination# = t * g1# + fngmst#(1900 + epochyear%(j%))
  657.   GreenwichDeclination# = (GreenwichDeclination# - INT(GreenwichDeclination#)) * pi2#
  658.   SinGD# = SIN(GreenwichDeclination#)
  659.   CosGD# = COS(GreenwichDeclination#)
  660.   ' Calculate sub-satellite position
  661.   SSPX# = (celestialX# * CosGD#) + (celestialY# * SinGD#)
  662.   SSPY# = -(celestialX# * SinGD#) + (celestialY# * CosGD#)
  663.   SSPZ# = celestialZ#
  664.   longitude! = 360 - fndegrees(arctan(SSPY#, SSPX#))
  665.   latitude! = fndegrees(fnarcsin(SSPZ# / r#))
  666.   XRange# = (SSPX# - GSxcoord)
  667.   YRange# = (SSPY# - GSycoord)
  668.   ZRange# = (SSPZ# - GSzcoord)
  669.   range# = SQR(XRange# * XRange# + YRange# * YRange# + ZRange# * ZRange#)
  670.   ' Translate to ground station coordinate system
  671.   GSCSZ# = (XRange# * CosGSLong# * CosGSLat#) + (YRange# * SinGSLong# * CosGSLat#) + (ZRange# * SinGSLat#)
  672.   GSCSX# = -(XRange# * CosGSLong# * SinGSLat#) - (YRange# * SinGSLong# * SinGSLat#) + (ZRange# * CosGSLat#)
  673.   GSCSY# = (YRange# * CosGSLong#) - (XRange# * SinGSLong#)
  674.   elevation! = fndegrees(fnarcsin(GSCSZ# / range#))
  675.   azimuth! = fndegrees(arctan(GSCSY#, GSCSX#))
  676.   ' Speed to/from ground station
  677.   IF OldT# <> t THEN Vr = ((range# - OldRange#) / (t - OldT#)) / 86400 ELSE Vr = -900000000
  678.   OldRange# = range#
  679.   OldT# = t
  680. RETURN
  681.  
  682. 9000 END
  683.  
  684. FUNCTION arctan (y, x)
  685. ' Modified arctangent function.  Returns a value between 0 and 2 PI.
  686. SHARED pi#, pi2#
  687.  
  688. IF x > 0 AND y >= 0 THEN
  689.   arctan = ATN(y / x)
  690. ELSEIF x > 0 AND y < 0 THEN
  691.   arctan = ATN(y / x) + pi2#
  692. ELSEIF x < 0 THEN
  693.   arctan = ATN(y / x) + pi#
  694. ELSEIF y >= 0 THEN
  695.   arctan = pi# / 2
  696. ELSE arctan = 1.5# * pi#
  697. END IF
  698. END FUNCTION
  699.  
  700. DEFSNG A-Z
  701. ' Converts a date string to the julian day of the year
  702. FUNCTION datetojulian% (d$)
  703.  
  704. day% = VAL(MID$(d$, 4, 2))
  705.  
  706. SELECT CASE LEFT$(d$, 2)
  707. CASE "01"
  708. CASE "02"
  709.   day% = day% + 31
  710. CASE "03"
  711.   day% = day% + 59
  712. CASE "04"
  713.   day% = day% + 90
  714. CASE "05"
  715.   day% = day% + 120
  716. CASE "06"
  717.   day% = day% + 151
  718. CASE "07"
  719.   day% = day% + 181
  720. CASE "08"
  721.   day% = day% + 212
  722. CASE "09"
  723.   day% = day% + 243
  724. CASE "10"
  725.   day% = day% + 273
  726. CASE "11"
  727.   day% = day% + 304
  728. CASE "12"
  729.   day% = day% + 334
  730. END SELECT
  731.  
  732. IF fnleap%(VAL(RIGHT$(d$, 4))) AND LEFT$(d$, 2) > "02" THEN
  733.   datetojulian% = day% + 1
  734. ELSE
  735.   datetojulian% = day%
  736. END IF
  737.  
  738. END FUNCTION
  739.  
  740. SUB getkey (k$)
  741.  
  742. k$ = ""
  743. WHILE k$ = ""
  744.         k$ = UCASE$(INKEY$)
  745. WEND
  746.  
  747. END SUB
  748.  
  749. ' Gets an input line, trims leading and trailing blanks
  750. ' then converts to upper-case before returning.
  751. FUNCTION gets$
  752. LINE INPUT i$
  753. gets$ = UCASE$(trim$(i$))
  754. END FUNCTION
  755.  
  756. ' Converts a julian day and year to a date string.
  757. FUNCTION juliantodate$ (julianday%, year%)
  758.  
  759. day% = INT(julianday%)
  760. IF day% = 60 THEN
  761.   IF fnleap%(year%) THEN tempdate$ = "02-29" ELSE tempdate$ = "03-01"
  762.   GOTO appendyear
  763. ELSE IF fnleap%(year%) AND day% > 60 THEN day% = day% - 1
  764. END IF
  765.  
  766. SELECT CASE day%
  767. CASE 1 TO 31
  768.   tempdate$ = "01-" + MID$(STR$(day% + 100), 3)
  769. CASE 32 TO 59
  770.   tempdate$ = "02-" + MID$(STR$(day% - 31 + 100), 3)
  771. CASE 60 TO 90
  772.   tempdate$ = "03-" + MID$(STR$(day% - 59 + 100), 3)
  773. CASE 91 TO 120
  774.   tempdate$ = "04-" + MID$(STR$(day% - 90 + 100), 3)
  775. CASE 121 TO 151
  776.   tempdate$ = "05-" + MID$(STR$(day% - 120 + 100), 3)
  777. CASE 152 TO 181
  778.   tempdate$ = "06-" + MID$(STR$(day% - 151 + 100), 3)
  779. CASE 182 TO 212
  780.   tempdate$ = "07-" + MID$(STR$(day% - 181 + 100), 3)
  781. CASE 213 TO 243
  782.   tempdate$ = "08-" + MID$(STR$(day% - 212 + 100), 3)
  783. CASE 244 TO 273
  784.   tempdate$ = "09-" + MID$(STR$(day% - 243 + 100), 3)
  785. CASE 274 TO 304
  786.   tempdate$ = "10-" + MID$(STR$(day% - 273 + 100), 3)
  787. CASE 305 TO 334
  788.   tempdate$ = "11-" + MID$(STR$(day% - 304 + 100), 3)
  789. CASE 335 TO 365
  790.   tempdate$ = "12-" + MID$(STR$(day% - 334 + 100), 3)
  791. END SELECT
  792.  
  793. appendyear: juliantodate$ = tempdate$ + "-" + MID$(STR$(year%), 2)
  794.  
  795. END FUNCTION
  796.  
  797. SUB plotloc (longitude!, latitude!, xcoord%, ycoord%)
  798. ' Returns the screen x and y coordinates for a given map location
  799.  
  800. ycoord% = CINT(.7111 * (90 - latitude!) + 3)
  801. IF longitude! >= 0 AND longitude! <= 270 THEN
  802.   xcoord% = CINT(477 - longitude! * 1.7444)
  803. ELSE xcoord% = CINT(1105 - longitude! * 1.7444)
  804. END IF
  805.  
  806. END SUB
  807.  
  808. SUB readgrnd
  809. SHARED gsid$, gsloc$, gslong!, gslat!, gsheight%, tzcor!
  810.  
  811. ' Format of USATGRND.DAT:
  812. '   gsid$    = Station Identifier
  813. '   gsloc$   = Station location (name)
  814. '   gslong!   = Longitude in degrees
  815. '   gslat!    = Latitude in degrees
  816. '   gsheight% = Height above sea level in meters
  817. '   tzcor!    = Time zone correction to UTC
  818. OPEN "USATGRND.DAT" FOR INPUT AS #1
  819. LINE INPUT #1, gsid$
  820. LINE INPUT #1, gsloc$
  821. INPUT #1, gslong!, gslat!, gsheight%, tzcor!
  822. CLOSE #1
  823.  
  824. END SUB
  825.  
  826. FUNCTION readtle% (tlefile$, name$(), epochyear%(), epochjulianday#(), inclination#(), epochRAAN#(), eccentricity#(), epochArgOfPerigee#(), epochMA#(), MeanMotion#(), epochrevolution&())
  827. ' ****** ESTABLISH SATELLITE ELEMENT MATRIX ******
  828. ' tlefile$ conforms to the NASA/NORAD three-line element set format.
  829.  
  830. OPEN tlefile$ FOR INPUT AS #1
  831.  
  832. FOR i% = 1 TO UBOUND(name$)
  833.   IF EOF(1) THEN EXIT FOR
  834.   LINE INPUT #1, name$(i%)
  835.   IF trim$(name$(i%)) = "" THEN EXIT FOR
  836.   LINE INPUT #1, y3$
  837.   LINE INPUT #1, t0$
  838.   epochyear%(i%) = VAL(MID$(y3$, 19, 2))
  839.   epochjulianday#(i%) = VAL(MID$(y3$, 21, 12))
  840.   inclination#(i%) = VAL(MID$(t0$, 9, 8))
  841.   epochRAAN#(i%) = VAL(MID$(t0$, 18, 8))
  842.   eccentricity#(i%) = VAL("." + MID$(t0$, 27, 7))
  843.   epochArgOfPerigee#(i%) = VAL(MID$(t0$, 35, 8))
  844.   epochMA#(i%) = VAL(MID$(t0$, 44, 8))
  845.   MeanMotion#(i%) = VAL(MID$(t0$, 53, 11))
  846.   epochrevolution&(i%) = VAL(MID$(t0$, 64, 5))
  847. NEXT
  848.  
  849. CLOSE #1
  850. readtle% = i% - 1
  851.            
  852. END FUNCTION
  853.  
  854. ' Converts an input date and time string from one time zone to another.
  855. SUB timezone (d0$, t0$, d1$, t1$, tzcor!)
  856.  
  857. julianday% = datetojulian(d0$)
  858.  
  859. year% = VAL(RIGHT$(d0$, 4))
  860.  
  861. hour! = VAL(LEFT$(t0$, 2)) + INT(tzcor!) + (VAL(MID$(t0$, 4, 2)) + tzcor! - INT(tzcor!)) / 60
  862. SELECT CASE hour!
  863. CASE IS < 0
  864.   hour! = hour! + 24
  865.   julianday% = julianday% - 1
  866. CASE IS >= 24
  867.   hour! = hour! - 24
  868.   julianday% = julianday% + 1
  869. END SELECT
  870. t1$ = fnt$(INT(hour!)) + ":" + fnt$(INT(60 * (hour! - INT(hour!)))) + MID$(t0$, 6)
  871.  
  872. SELECT CASE julianday%
  873. CASE 0
  874.   year% = year% - 1
  875.   IF fnleap%(year%) THEN julianday% = 366 ELSE julianday% = 365
  876. CASE IS > 365
  877.   IF NOT fnleap%(year%) THEN julianday% = 1: year% = year% + 1
  878. END SELECT
  879.  
  880. d1$ = juliantodate$(julianday%, year%)
  881. END SUB
  882.  
  883. ' Trims leading and trailing spaces from a string.
  884. FUNCTION trim$ (i$)
  885.           
  886. trim$ = RTRIM$(LTRIM$(i$))
  887.  
  888. END FUNCTION
  889.  
  890.